﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace NCM_MSTest.Alg
{
    public delegate double Fxy(double x, double y);
    public delegate double Fxyz(double x, double y, double z);

    /// <summary>
    /// 常微分方程的数值解
    /// 常微分方程的定解问题分为：初值问题和边值问题
    ///     求解初值问题：求同时满足方程和初始条件的解
    ///     求解边值问题：求同时满足方程和边界条件的解
    /// 一阶常微分方程初值问题的三种数值解法：
    ///     欧拉方法和改进的欧拉方法、龙格库塔方法、阿达姆斯方法
    /// 边值问题的解法：
    ///     1.将边值问题转化为初值问题求解
    ///     2.利用数值微分公式，将常微分方程边值问题转化为差分方程
    /// </summary>
    internal class OrdinaryDifferentialEquation
    {
        /// <summary>
        /// 改进的欧拉方法
        /// 欧拉方法：
        ///     y_{i+1} = y_i + h*f(x_i,y_i),i=0,1,2...
        ///     具有一阶精度O（h^2）
        /// 改进的欧拉方法
        ///     y_{i+1} = y_i +h/2*[f(x_i,y_i) + f(x_{i+1},y_i+hf(x_i,y_i))],i=0,1,2...
        ///     具有二阶精度O（h^3）
        /// </summary>
        /// <param name="fxy"></param>
        /// <param name="iy"></param>
        /// <param name="ix"></param>
        /// <param name="h"></param>
        /// <param name="xs"></param>
        /// <param name="ys"></param>
        public static void ImprovedEulerMethod(Fxy fxy
            , double iy, double ix
            , double h
            , double[] xs
            , out double[] ys
            )
        {
            ys = new double[xs.Length];
            double xf = xs.Last();
            double x0 = ix;
            double y0 = iy;
            int n = (int)Math.Floor((xf - xs.First()) / h);
            double x = x0;
            double y = y0;
            double y1 = y;
            double eps1 = 0, eps2 = 0;
            Console.Write("{i}\t{x1}\t{y1}\t{y21}\t{y2}\t");
            Console.WriteLine("{ytrue}\t{eps1%}\t{eps2%}");
            for (int i = 0; i <= n; i++)
            {
                y1 = y1 + h * fxy(x, y1);//欧拉法
                double y21 = y + h * fxy(x, y);
                double x1 = x + h;
                double y2 = y + 0.5 * h * (fxy(x, y) + fxy(x1, y21));//改进的欧拉法
                double ytrue = Math.Pow(Math.E, -x1) + x1;
                eps1 = Math.Abs((ytrue - y1) / ytrue) * 100.0;
                eps2 = Math.Abs((ytrue - y2) / ytrue) * 100.0;
                Console.Write($"{i + 1}\t{x1:f5}\t{y1:f5}\t{y21:f5}\t{y2:f5}\t");
                Console.WriteLine($"{ytrue:f5}\t{eps1:f5}\t{eps2:f5}");
                x = x1;
                y = y2;
                ys[i] = y2;
            }
        }

        /// <summary>
        /// 龙格库塔法
        /// 在积分小区域内多预报几个点的斜率，然后将其线性组合作为平均斜率的近似值
        /// 二阶龙格库塔公式
        ///     x_{i+p}=x_i+ph
        ///     y_{i+1}=y_i+h*(L1k1+L2k2)
        ///     k1=f(xi,yi)
        ///     y_{i+p}=y_i+p*h*k1
        ///     k2=f(x_{i+p},y_{i+p})
        /// 三阶龙格库塔公式
        ///     ...
        /// 四阶龙格库塔公式
        /// </summary>
        /// <param name="f"></param>
        public static void Runge_Kutta(Fxy f
            , double iy, double ix
            , double h
            , double[] xs
            , out double[] ys
            )
        {
            ys = new double[xs.Length];
            double xf = xs.Last();
            double x0 = ix;
            double y0 = iy;
            double x = x0;
            double y = y0;
            double fk1, fk2, fk3, fk4;
            double ytrue;
            double eps;
            Console.WriteLine("{x:f6}\t {y:f6}\t {ytrue:f6}\t {eps:e6}");
            while (x - xf < 0)
            {
                fk1 = f(x, y);
                fk2 = f(x + h / 2, y + fk1 * h / 2);
                fk3 = f(x + h / 2, y + fk2 * h / 2);
                fk4 = f(x + h, y + fk3 * h);
                y = y + (fk1 + 2.0 * (fk2 + fk3) + fk4) * h / 6.0;
                x = x + h;
                ytrue = Math.Pow(Math.E, -x) + x;
                eps = Math.Abs((ytrue - y) / ytrue) * 100.0;
                Console.WriteLine($"{x:f6}\t {y:f6}\t {ytrue:f6}\t {eps:e6}");
            }
        }

        /// <summary>
        /// 阿达姆斯方法 预估-校正方法
        /// </summary>
        /// <param name="f"></param>
        /// <param name="iy"></param>
        /// <param name="ix"></param>
        /// <param name="h"></param>
        /// <param name="xs"></param>
        /// <param name="ys"></param>
        public static void Adams(Fxy f
            , double iy, double ix
            , double h
            , double[] xs
            , out double[] ys
            )
        {
            ys = new double[xs.Length];
            double xf = xs.Last();
            double x0 = ix;
            double y0 = iy;
            double fk1, fk2, fk3, fk4;
            double ytrue;
            double eps;
            int n = (int)Math.Round((xf - x0) / h);
            double[] x = new double[n + 1];
            double[] y = new double[n + 1];
            for (int i = 0; i <= n; i++)
            {
                x[i] = x0 + i * h;
            }
            y[0] = y0;
            Console.WriteLine("{x[i]:f6}\t {y[i]:f6}\t {ytrue:f6}\t {eps:e6}");
            Console.WriteLine("starting values from Runge-kutta:");
            for (int i = 1; i < 4; i++)
            {
                fk1 = f(x[i - 1], y[i - 1]);
                fk2 = f(x[i - 1] + h / 2, y[i - 1] + fk1 * h / 2);
                fk3 = f(x[i - 1] + h / 2, y[i - 1] + fk2 * h / 2);
                fk4 = f(x[i - 1] + h, y[i - 1] + fk3 * h);
                y[i] = y[i - 1] + (fk1 + 2.0 * (fk2 + fk3) + fk4) * h / 6;
                ytrue = Math.Pow(Math.E, -x[i]) + x[i];
                eps = Math.Abs((ytrue - y[i]) / ytrue) * 100.0;
                Console.WriteLine($"{x[i]:f6}\t {y[i]:f6}\t {ytrue:f6}\t {eps:e6}");
            }
            Console.WriteLine("continuation of solution by admas method:");
            for (int i = 3; i < n; i++)
            {
                y[i + 1] = y[i] + h * (55 * f(x[i], y[i]) - 59 * f(x[i - 1], y[i - 1])
                    + 37 * f(x[i - 2], y[i - 2])
                    - 9 * f(x[i - 3], y[i - 3])) / 24;
                y[i + 1] = y[i] + h * (9 * f(x[i + 1], y[i + 1]) + 19 * f(x[i], y[i])
                    - 5 * f(x[i - 1], y[i - 1]) + f(x[i - 2], y[i - 2])) / 24;
                ytrue = Math.Pow(Math.E, -x[i + 1]) + x[i + 1];
                eps = Math.Abs((ytrue - y[i + 1]) / ytrue) * 100.0;
                Console.WriteLine($"{x[i + 1]:f6}\t {y[i + 1]:f6}\t {ytrue:f6}\t {eps:e6}");
            }
        }

        /// <summary>
        /// 求解一阶常微分方程组
        /// 四阶龙格-库塔方法
        /// </summary>
        /// <param name="f"></param>
        /// <param name="g"></param>
        /// <param name="x0"></param>
        /// <param name="y10"></param>
        /// <param name="y20"></param>
        /// <param name="xf"></param>
        /// <param name="h"></param>
        /// <param name="ys"></param>
        public static void Runge_kutta_4(Fxyz f, Fxyz g
            , double x0
            , double y10
            , double y20
            , double xf
            , double h
            , double[] ys
            )
        {
            double x = x0;
            double y1 = y10;
            double y2 = y20;
            double fk1, fk2, fk3, fk4;
            double fl1, fl2, fl3, fl4;
            Console.WriteLine("{x:f6}\t {y1:f6}\t {y2:f6}\t {x - xf}");
            while (Math.Abs(x - xf) > 1e-15)
            {
                fk1 = f(x, y1, y2);
                fl1 = g(x, y1, y2);
                fk2 = f(x + h / 2, y1 + fk1 * h / 2, y2 + fl1 * h / 2);
                fl2 = g(x + h / 2, y1 + fk1 * h / 2, y2 + fl1 * h / 2);
                fk3 = f(x + h / 2, y1 + fk2 * h / 2, y2 + fl2 * h / 2);
                fl3 = g(x + h / 2, y1 + fk2 * h / 2, y2 + fl2 * h / 2);
                fk4 = f(x + h, y1 + fk3 * h, y2 + fl3 * h);
                fl4 = g(x + h, y1 + fk3 * h, y2 + fl3 * h);

                y1 = y1 + (fk1 + 2.0 * (fk2 + fk3) + fk4) * h / 6;
                y2 = y2 + (fl1 + 2.0 * (fl2 + fl3) + fl4) * h / 6;
                x = x + h;
                Console.WriteLine($"{x:f6}\t {y1:f6}\t {y2:f6}\t {x - xf}");
            }
        }

        /// <summary>
        /// 二阶常微分方程四阶龙格库塔公式计算
        /// </summary>
        public static void SecondOrder_Runge_kutta_4(Fxy f
            , double t0, double tf
            , double y0, double v0
            , double h
            )
        {
            double t = t0;
            double y = y0;
            double v = v0;
            double a = f(t, v);
            int n0 = 0;
            double fl1 = 0;
            double fl2 = 0;
            double fl3 = 0;
            double fl4 = 0;
            Console.WriteLine("{n0}\t {t:f6}\t {y:f6}\t  {v:f6}\t {a:f6}");
            while (t <= tf)
            {
                fl1 = f(t, v);
                fl2 = f(t + h / 2, v + fl1 * h / 2);
                fl3 = f(t + h / 2, v + fl2 * h / 2);
                fl4 = f(t + h, v + fl3 * h);
                y = y + h * v + (fl1 + fl2 + fl3) * h * h / 6.0;
                v = v + (fl1 + 2.0 * (fl2 + fl3) + fl4) * h / 6.0;
                t = t + h;
                a = f(t, v);
                n0 = n0 + 1;
                if (n0 % 200 == 0)
                    Console.WriteLine($"{n0}\t {t:f6}\t {y:f6}\t  {v:f6}\t {a:f6}");
            }
        }

        /// <summary>
        /// 常微分方程边值问题的解法
        /// 二阶线性常微分方程：y''+p(x)y'+q(x)y=f(x)
        /// 第一边值条件：y(a)=alpha,y(b)=beta
        /// 第二边值条件：y'(a)=alpha,y'(b)=beta
        /// 第三边值条件：y'(a) - b_1*y(a)=alpha,y'(b)-b_2*y(b)=beta
        /// 本方法计算第一边界条件
        /// 计算过程：
        ///     1.将边值问题转化为初值问题
        ///         求解初值问题：y''+p(x)y'+q(x)y=f(x),(x>a), y(a)=alpha,y'(a)=0
        ///             其数值解为：y^1(x_i)
        ///         求解初值问题：y''+p(x)y'+q(x)y=0,(x>a),y(a)=0,y'(a)=1
        ///             其数值解为：y^2(x_i)
        ///     2.则原边值问题的通解：y(x_i)=y^1(x_i) + Cy^2(x_i)
        ///         y(b)=y^1(x_n)+Cy^2(x_n)=Beta，确定常数C，即得到边值问题的同解
        ///     F1 = f(x) - [p(x)y'+q(x)y]
        ///     F2 = 0    - [p(x)y'+q(x)y]
        /// </summary>
        public static void Runge_kutta_Border(Fxyz f1, Fxyz f2
            , double h
            , double x0, double xf
            , double y110
            , double y120
            , double y210
            , double y220
            , double y20
            , out double[] y
            )
        {
            int n = (int)Math.Floor((xf - x0) / h);
            y = new double[n];
            int i = 0;
            double x = x0;
            double y11 = y110, y12 = y120;
            double y21 = y210, y22 = y220;
            double[] y1 = new double[n];
            double[] y2 = new double[n];
            Console.WriteLine($"N:{n}");
            Console.WriteLine($"{x:f2}\t{y11:f5}\t{y12:f5}\t{y21:f5}\t{y22:f5}");
            while (x < xf && i < n)
            {
                Rkutta(f1, x, h, ref y11, ref y12);
                Rkutta(f2, x, h, ref y21, ref y22);
                y1[i] = y11;
                y2[i] = y21;
                x += h;
                i++;
                Console.WriteLine($"{i} {x:f2}\t{y11:f5}\t{y12:f5}\t{y21:f5}\t{y22:f5}");
            }
            double C = (y20 - y1[n - 1]) / y2[n - 1];
            Console.WriteLine($"C:{C}");
            Console.WriteLine("The results is:");
            for (i = 0; i < n; i++)
            {
                y[i] = y1[i] + C * y2[i];
                Console.WriteLine($"{i + 1} \t {x0 + (i + 1) * h:f2} \t {y[i]}");
            }
        }

        private static void Rkutta(Fxyz f
            , double x, double h
            , ref double y1
            , ref double y2
            )
        {
            double fl1 = f(x, y1, y2);
            double fl2 = f(x + h / 2, y1 + y2 * h / 2, y2 + fl1 * h / 2);
            double fl3 = f(x + h / 2, y1 + y2 * h / 2 + fl1 * h * h / 4, y2 + fl2 * h / 2);
            double fl4 = f(x + h, y1 + y2 * h + fl2 * h * h / 2, y2 + fl3 * h);

            y1 = y1 + h * y2 + (fl1 + fl2 + fl3) * h * h / 6;
            y2 = y2 + (fl1 + 2 * (fl2 + fl3) + fl4) * h / 6;
        }

        /// <summary>
        /// 边值问题的差分解法
        /// y''+p(x)y'+q(x)y=f(x)
        /// y(a)=alpha,y(b)=beta
        /// 
        /// </summary>
        public static void DifferenceMethod_TwoOrderBoundary(int n, double h, double ai, double bf
            , double[] a, double[] b, double[] c, double[] d)
        {
            for (int i = 0; i < n; i++)
            {
                int ti = i;
                d[i] = 2.0 * h * h * (-ti * h) - 4;
                if (i < n)
                {
                    a[i] = 2.0 - h * h * (ti + 1);
                    c[i] = 2.0 + h * h * ti;
                    if (i == 0)
                    {
                        b[i] = 4.0 * h * h * h - (2.0 - h * h) * ai;
                    }
                    else
                    {
                        b[i] = 4.0 * h * h * h * ti;
                    }
                }
                else
                {
                    b[i] = 4.0 * h * h * h * ti - (2.0 + h * h * ti) * bf;
                }
                Trid(n, a, d, c, b);
                Console.WriteLine($"{i},{i * h},{b[i]}");
            }
        }

        private static void Trid(int n, double[] a, double[] d, double[] c, double[] b)
        {
            double[] a1 = new double[n];
            double[] b1 = new double[n];
            for (int i = 1; i < n; i++)
            {
                int m = i - 1;
                a1[m] = c[m] / d[m];
                b1[m] = b[m] / d[m];
                d[i] = d[i] - a1[m] * a[m];
                b[i] = b[i] - b1[m] * a[m];
            }
            b[n - 1] = b[n - 1] / d[n - 1];
            for (int i = 0; i < n - 2; i++)
            {
                int k = n - 2 - i;
                b[k] = b1[k] - a1[k] * b[k + 1];
            }
        }
    }
}
